//2018-4-13
#ifndef WGT_H
#define WGT_H
#include <vector>
#include <string>
#include <iostream>
#include "gdal_priv.h"
#include "ogr_spatialref.h"
#include <cmath>
#include <cstring>
#include <fstream>
#include <math.h>

using std::cout ;
using std::vector ;
using std::string ;
using std::endl ;

//判断文件是否存在
bool test_file_exists(const std::string& name) {
  std::ifstream f(name.c_str());
  return f.good();
}

struct ccPoint
{
	float x, y;
};

struct ccRasterInfo
{
	double trans[6];
	int xsize, ysize;
	std::string proj4str;
	std::string projstr;
};

struct ccDataInfo
{
	double valid0 , valid1 ;
	double scale , offset ;
	inline ccDataInfo(){} ;
	inline ccDataInfo(double v0,double v1,double sc,double off):valid0(v0),valid1(v1),scale(sc),offset(off){};
};

class wgt_CoorConverter
{
	std::vector<std::string> m_predefinedProj4;
	std::vector<std::string> m_predefinedCode;

	OGRCoordinateTransformation* m_pConverter;
	OGRCoordinateTransformation* m_pRConverter;
public:
	std::string m_proj4str;
	wgt_CoorConverter();
	~wgt_CoorConverter();
	bool prepare(std::string proj4str);
	void convert(double& xlon, double& ylat);
	///convert other coordinate into longlat
	bool rconvert(double& x,double& y );
	
	void preparePredefined(std::string code);
	void printPredefined();

	void makeLongLatGridAndLabels(std::string gridfile, std::string labelfile , 
		double lon0 , double lon1 , double lat0 , double lat1 );
	void convertTextfile(std::string infile, std::string outfile);
	int splitline(std::string line, double& v0, double& v1, std::string& str);
	void getCoorRange(std::string file, double& x0, double& x1, double& y0, double& y1);
	bool getPixelRangeByLongLat(std::string file,
		double lon0, double lon1, double lat0, double lat1,
		int& x0, int& x1, int& y0, int& y1);
	std::string getProj4(GDALDataset* ds);
	bool datacombineAve(std::vector<std::string>& infiles, std::string outfile, double valid0, double valid1, std::string& error, bool ignorebadfile);
	bool datacombineMax(std::vector<std::string>& infiles, std::string outfile, double valid0, double valid1, std::string& error, bool ignorebadfile);
	bool datacombineSum(std::vector<std::string>& infiles, std::string outfile, std::string& error, bool ignorebadfile);
	bool datacombineAve2(std::vector<std::string>& infiles, std::string outfile, double valid0, double valid1, std::string& error, bool ignorebadfile);
	bool datacombineAve3(std::vector<std::string>& infiles, std::string outfile, double valid0, double valid1, std::string& error, bool ignorebadfile, vector<bool>& hasvcntVec);
	ccRasterInfo getRasterInfo(std::string filepath);
	bool dataSubtract(std::string fileA , std::string fileB , std::string outfile , 
			ccDataInfo ia, ccDataInfo ib, double outfillvalue, std::string& error ) ;
};

wgt_CoorConverter::wgt_CoorConverter()
{
	m_pConverter = 0;
	m_pRConverter = 0;

	m_predefinedCode.push_back("geos1047");
	m_predefinedProj4.push_back("+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=104.7 +no_defs");

	m_predefinedCode.push_back("geos0995");
	m_predefinedProj4.push_back("+proj=geos +h=35785863 +a=6378137.0 +b=6356752.3 +lon_0=99.5 +no_defs");

	m_predefinedCode.push_back("npolar095");
	m_predefinedProj4.push_back("+proj=stere +lat_0=90 +lat_ts=0 +lon_0=95 +k=1 +x_0=0 +y_0=0 +a=6371200 +b=6371200 +units=m +no_defs");

	m_predefinedCode.push_back("albers105");
	m_predefinedProj4.push_back("+proj=aea +ellps=krass +lon_0=105 +lat_1=25 +lat_2=47");
}

wgt_CoorConverter::~wgt_CoorConverter()
{
	if (m_pConverter)
	{
		OCTDestroyCoordinateTransformation(m_pConverter);
		m_pConverter = 0;
	}
	if (m_pRConverter)
	{
		OCTDestroyCoordinateTransformation(m_pRConverter);
		m_pConverter = 0;
	}
}

bool wgt_CoorConverter::prepare(std::string proj4str)
{
	if (m_pConverter)
	{
		OCTDestroyCoordinateTransformation(m_pConverter);
		m_pConverter = 0;
	}
	if (m_pRConverter)
	{
		OCTDestroyCoordinateTransformation(m_pRConverter);
		m_pConverter = 0;
	}
	m_proj4str = proj4str;
	OGRSpatialReference wgs84;
	wgs84.SetWellKnownGeogCS("WGS84");
	char proj2[2048];
	sprintf(proj2, "%s", m_proj4str.c_str());
	OGRSpatialReference targetSr;
	targetSr.importFromProj4(proj2);
	m_pConverter = OGRCreateCoordinateTransformation(&wgs84,
		&targetSr);
	m_pRConverter = OGRCreateCoordinateTransformation(&targetSr,
		&wgs84);
	return true;
}

///lon,lat-->x,y
void wgt_CoorConverter::convert(double& xlon, double& ylat)
{
	if (m_pConverter)
	{
		m_pConverter->Transform(1, &xlon, &ylat);
	}
}
///x,y-->lon,lat
bool wgt_CoorConverter::rconvert(double& x, double& y)
{
	if (m_pRConverter)
	{
		bool isok = m_pRConverter->Transform(1, &x, &y);
		if (isok)
		{
			return true;
		}
		else
		{
			return false;
		}
	}
	else
	{
		return false;
	}
}

void wgt_CoorConverter::preparePredefined(std::string code)
{
	bool ok = false;
	for (int i = 0; i < m_predefinedCode.size(); ++i)
	{
		if (code == m_predefinedCode[i])
		{
			ok = true;
			this->prepare(m_predefinedProj4[i]);
			break;
		}
	}
	if (ok == false)
	{
		std::cout << "wgt_CoorConverter can not find predefined proj4 string : " << code << std::endl;
	}
}
void wgt_CoorConverter::printPredefined()
{
	for (int i = 0; i < m_predefinedCode.size(); ++i)
	{
		std::cout << i << "\t" << m_predefinedCode[i] << "\t" << m_predefinedProj4[i] << std::endl;
	}
}

void wgt_CoorConverter::makeLongLatGridAndLabels(std::string gridfile, std::string labelfile, double lon0, double lon1, double lat0, double lat1)
{
	lon0 = floor(lon0+0.5);
	lon1 = floor(lon1+0.5);
	lat0 = floor(lat0+0.5);
	lat1 = floor(lat1+0.5);

	double dlon = lon1 - lon0;
	if (dlon < 5)
	{
		dlon = 1;
	}
	else if (dlon < 20)
	{
		dlon = 5;
	}
	else if (dlon < 50 )
	{
		dlon = 10;
	}
	else
	{
		dlon = 20;
	}

	double dlat = lat1 - lat0;
	if (dlat < 5)
	{
		dlat = 1;
	}
	else if (dlat < 20)
	{
		dlat = 5;
	}
	else if (dlat < 50)
	{
		dlat = 10;
	}
	else
	{
		dlat = 20;
	}

	FILE* pf = fopen(gridfile.c_str(), "w");
	FILE* pf1 = fopen(labelfile.c_str(), "w");
	fprintf(pf, "#x  y\n");
	fprintf(pf1, "#x  y  label\n");
	for (double tlon = lon0; tlon <= lon1; tlon += dlon)
	{
		bool labeled = false;
		for (double tlat = lat0 ; tlat <= lat1; tlat += 0.05)
		{
			if (!labeled)
			{
				labeled = true;
				double labellon = tlon + dlon;
				int tloni = tlon ;
				int dloni = dlon ;
				if (labellon < lon1)
				{
					double labellat = tlat + dlat / 2;
					double x1 = labellon;
					double y1 = labellat;
					this->convert(x1, y1);
					if (labellon > 0)
					{
						fprintf(pf1, "%f %f %8.1f\\260E\n", x1, y1, labellon);
					}
					else if (labellon == 0)
					{
						fprintf(pf1, "%f %f %8.1f\\260\n", x1, y1, labellon);
					}
					else
					{
						fprintf(pf1, "%f %f %8.1f\\260W\n", x1, y1, labellon);
					}
				}
				
			}
			double x = tlon;
			double y = tlat;
			this->convert(x, y);
			fprintf(pf, "%f  %f\n", x, y);
		}
		fprintf(pf, "\n");
	}
	for (double tlat = lat0; tlat <= lat1; tlat += dlat )
	{
		bool labeled = false;
		for (double tlon = lon0; tlon <= lon1; tlon += 0.05)
		{
			if (!labeled)
			{
				labeled = true;
				double labellat = tlat ;
				if (labellat < lat1)
				{
					double labellon = tlon + dlon / 2;
					double x1 = labellon;
					double y1 = labellat;
					this->convert(x1, y1);
					if (labellat > 0)
					{
						fprintf(pf1, "%f %f %8.1f\\260N\n", x1, y1, labellat);
					}
					else if (labellat == 0)
					{
						fprintf(pf1, "%f %f %8.1f\\260\n", x1, y1, labellat);
					}
					else
					{
						fprintf(pf1, "%f %f %8.1f\\260S\n", x1, y1, labellat);
					}
				}

			}
			double x = tlon;
			double y = tlat;
			this->convert(x, y);
			fprintf(pf, "%f  %f\n", x, y);
		}
		fprintf(pf, "\n");
	}
	fclose(pf);
	fclose(pf1);
}



//2018-1-17 一行前两个数字保留，剩下的作为整个字符串
int wgt_CoorConverter::splitline(std::string line, double& v0, double& v1,std::string& str)
{
	str = "";
	int len = line.length();
	std::string part;
	int nreal = 0;
	double vtwo[2] = { 0.0,0.0 };
	int istr = len ;
	for (int i = 0; i < len; ++i)
	{
		if (nreal < 2)
		{
			bool isdig = false;
			if (line[i] == '-' || line[i] == '.' || (line[i] >= '0' && line[i] <= '9'))
			{
				part.push_back(line[i]);
				isdig = true;
			}
			
			if( (i==len-1 || isdig==false ) && part.length()>0 )
			{
				double tval = atof(part.c_str());
				part = "";
				vtwo[nreal] = tval;
				++nreal;
			}
		}
		else
		{
			istr = i;
			break;
		}
	}
	v0 = vtwo[0];
	v1 = vtwo[1];
	if (istr < len)
	{
		str = line.substr(istr);
	}
	

	return nreal ;
}

//2018-1-17
void wgt_CoorConverter::convertTextfile(std::string infile, std::string outfile)
{
	FILE* pf1 = fopen(outfile.c_str(), "w");
	std::ifstream infs(infile.c_str());
	std::string line;
	double v0, v1;
	std::string vstr;
	while (std::getline(infs, line))
	{
		int nreal = this->splitline(line, v0, v1, vstr);
		if (nreal == 2)
		{
			this->convert(v0, v1);
			fprintf(pf1, "%f %f %s\n", v0, v1, vstr.c_str());
		}
		else
		{
			fprintf(pf1, "%s\n", line.c_str());
		}
	}
	infs.close();
	fclose(pf1);
}

//2017-1-17
void wgt_CoorConverter::getCoorRange(std::string file, double& x0, double& x1, double& y0, double& y1)
{
	GDALDataset* ds = (GDALDataset*) GDALOpen(file.c_str(), GA_ReadOnly);
	if (ds)
	{
		double trans[6];
		ds->GetGeoTransform(trans);
		x0 = trans[0];
		y0 = trans[3] + ds->GetRasterYSize() * trans[5];
		x1 = trans[0] + ds->GetRasterXSize() * trans[1];
		y1 = trans[3];

		GDALClose(ds);
	}
}

//2017-1-17
std::string wgt_CoorConverter::getProj4( GDALDataset* ds )
{
	OGRSpatialReference osr;
	char prjbuff[2048];
	char* ptrPrj = prjbuff;
	char* ptr2 = 0;
	strcpy(prjbuff, ds->GetProjectionRef());
	osr.importFromWkt(&ptrPrj);
	osr.exportToProj4(&ptr2);
	std::string proj4str = std::string(ptr2);
	CPLFree(ptr2); ptr2 = 0;
	return proj4str;
}


//2017-1-17
bool wgt_CoorConverter::getPixelRangeByLongLat(std::string file,
	double lon0,double lon1,double lat0,double lat1,
	int& x0, int& x1, int& y0, int& y1)
{
	GDALDataset* ds = (GDALDataset*)GDALOpen(file.c_str(), GA_ReadOnly);
	if (ds)
	{
		double trans[6];
		ds->GetGeoTransform(trans);
		int xsize = ds->GetRasterXSize();
		int ysize = ds->GetRasterYSize();
		std::string prj4 = this->getProj4(ds);
		GDALClose(ds); ds = 0;
		this->prepare(prj4);
		x0 = 99999;
		x1 = 0;
		y0 = 99999;
		y1 = 0;
		double latarr[2] = { lat0,lat1 };
		for (double tlon = lon0; tlon < lon1 + 1; tlon += 1)
		{
			for (int j = 0; j < 2; ++j)
			{
				double tx1 = tlon;
				double ty1 = latarr[j];
				this->convert(tx1, ty1);
				int pixelx = floor( (tx1 - trans[0]) / trans[1] +0.5) ;
				int pixely = floor((ty1 - trans[3]) / trans[5] +0.5);
				x0 = MIN(pixelx, x0);
				x1 = MAX(pixelx, x1);
				y0 = MIN(pixely, y0);
				y1 = MAX(pixely, y1);
			}
		}
		x0 = MAX(0, x0);
		x1 = MIN(xsize - 1, x1);
		y0 = MAX(0, y0);
		y1 = MIN(ysize - 1, y1);
		return true;
	}
	else
	{
		return false;
	}
}

//2018-1-22
bool wgt_CoorConverter::datacombineAve(std::vector<std::string>& infiles, 
	std::string outfile, double valid0, double valid1, std::string& error , bool ignorebadfile)
{
	cout<<"wgt_CoorConverter::datacombineAve v1.0"<<endl ;
	std::cout << "starting combine for " << outfile <<" ignorebadfile:"<<ignorebadfile << std::endl;
	//get output parameters.
	int xsize(0), ysize(0);
	GDALDataType datatype;
	double trans[6];
	char prjbuff[2048];
	for (int i = 0; i < infiles.size(); ++i)
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen(infiles[i].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			xsize = tds->GetRasterXSize();
			ysize = tds->GetRasterYSize();
			datatype = tds->GetRasterBand(1)->GetRasterDataType();
			strcpy(prjbuff, tds->GetProjectionRef());
			tds->GetGeoTransform(trans);
			GDALClose(tds);
			break;
		}
	}

	if (xsize == 0)
	{
		error = "Failed to get raster info,proj,geotrans.";
		return false;
	}

	int asize = xsize * ysize;
	float* dataBuffer = new float[asize];
	short* countBuffer = new short[asize];
	float* outBuffer = new float[asize];
	if (dataBuffer == 0
		|| countBuffer == 0
		|| outBuffer == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to alloc memories.";
		return false;
	}

	memset(countBuffer, 0, asize * sizeof(short));
	for (int ifile = 0; ifile < infiles.size() ; ++ifile )
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen( infiles[ifile].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			int nx = tds->GetRasterXSize();
			int ny = tds->GetRasterYSize();
			if (nx != xsize || ny != ysize) {
				std::cout << "*** Error : " << infiles[ifile] << " raster data size is not equal with others. out." << std::endl;
				if( ignorebadfile == false )
				{
					if (dataBuffer) delete[] dataBuffer;
					if (outBuffer) delete[] outBuffer;
					if (countBuffer) delete[] countBuffer;
					GDALClose(tds);
					error = "difference raster size " + infiles[ifile] ;
					return false ;
				}
			}
			else
			{
				tds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
					dataBuffer, xsize, ysize, GDT_Float32, 0, 0, 0);
				for (int it = 0; it < asize; ++it)
				{
					if (countBuffer[it] == 0)
					{
						outBuffer[it] = dataBuffer[it];
						if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
						{
							countBuffer[it]+=1;
						}
					}
					else
					{
						if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
						{
							outBuffer[it] += dataBuffer[it];
							countBuffer[it] += 1;
						}
					}
				}
			}
			GDALClose(tds);
		}
		else
		{
			std::cout << "Failed to open " + infiles[ifile] << std::endl;
			if( ignorebadfile == false )
			{
				if (dataBuffer) delete[] dataBuffer;
				if (outBuffer) delete[] outBuffer;
				if (countBuffer) delete[] countBuffer;
				error = "failed to open " + infiles[ifile] ;
				return false ;
			}
		}
		std::cout << ifile << "/" << infiles.size() << std::endl;
	}

	for (int it = 0; it < asize; ++it)
	{
		if (countBuffer[it] > 0)
		{
			outBuffer[it] = outBuffer[it] / countBuffer[it];
		}
	}

	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (driver == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to get driver of tiff.";
		return false;
	}
	GDALDataset* outds = driver->Create(outfile.c_str(), xsize, ysize, 1, datatype, 0);
	if (outds == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to create out put : "+outfile ;
		return false;
	}

	outds->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,outBuffer , xsize,ysize,GDT_Float32 , 0,0,0);
	{
		outds->SetProjection(prjbuff) ;
		outds->SetGeoTransform(trans) ;
	}
	
	delete[] dataBuffer; dataBuffer = 0;
	delete[] countBuffer; countBuffer = 0;
	delete[] outBuffer; outBuffer = 0;

	GDALClose(outds);
	std::cout << "combine ok." << std::endl;
	return true;
}


//2018-4-11 有效个数得以保存 .vcnt.tif 计算输入文件需要全部都有.vcnt.tif或者全没有
bool wgt_CoorConverter::datacombineAve2(std::vector<std::string>& infiles, std::string outfile, double valid0, double valid1, std::string& error, bool ignorebadfile)
{
	cout<<"wgt_CoorConverter::datacombineAve v2.0"<<endl ;
	std::cout << "starting combine for " << outfile <<" ignorebadfile:"<<ignorebadfile << std::endl;

	vector<string> cntfiles ;
	for( int i = 0 ; i<infiles.size() ; ++ i )
	{
		string cfile1 = infiles[i] + ".vcnt.tif" ;
		if( test_file_exists(cfile1) )
		{
			cntfiles.push_back(cfile1) ;
		}
	}

	bool useVcnt = false ;
	if( cntfiles.size() == 0 )
	{
		useVcnt = false ;
		cout<<"not using vcnt."<<endl ;
	}else 
	{
		if( cntfiles.size() == infiles.size() )
		{
			useVcnt = true ;
			cout<<"using vcnt."<<endl ;
		}else
		{
			error = "Error: not all input files have a .vcnt.tif file.";
			return false;
		}
	}


	//get output parameters.
	int xsize(0), ysize(0);
	GDALDataType datatype;
	double trans[6];
	char prjbuff[2048];
	for (int i = 0; i < infiles.size(); ++i)
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen(infiles[i].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			xsize = tds->GetRasterXSize();
			ysize = tds->GetRasterYSize();
			datatype = tds->GetRasterBand(1)->GetRasterDataType();
			strcpy(prjbuff, tds->GetProjectionRef());
			tds->GetGeoTransform(trans);
			GDALClose(tds);
			break;
		}
	}

	if (xsize == 0)
	{
		error = "Failed to get raster info,proj,geotrans.";
		return false;
	}

	int asize = xsize * ysize;
	float* dataBuffer = new float[asize];
	short* countBuffer = new short[asize];
	float* outBuffer = new float[asize];
	short* vcntbuffer = 0 ;
	if( useVcnt )
	{
		vcntbuffer = new short[asize];
		if( vcntbuffer==0 )
		{
			error = "Failed to alloc memories for vcntbuffer.";
			if (dataBuffer) delete[] dataBuffer;
			if (outBuffer) delete[] outBuffer;
			if (countBuffer) delete[] countBuffer;
			
			return false;
		}
	}

	if (dataBuffer == 0
		|| countBuffer == 0
		|| outBuffer == 0 )
	{
		error = "Failed to alloc memories.";
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		if (vcntbuffer) delete[] vcntbuffer;
		
		return false;
	}

	memset(countBuffer, 0, asize * sizeof(short));
	for (int ifile = 0; ifile < infiles.size() ; ++ifile )
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen( infiles[ifile].c_str(), GA_ReadOnly);

		if (tds != 0 )
		{
			int nx = tds->GetRasterXSize();
			int ny = tds->GetRasterYSize();
			if (nx != xsize || ny != ysize) {
				std::cout << "*** Error : " << infiles[ifile] << " raster data size is not equal with others. out." << std::endl;
				if( ignorebadfile == false )
				{
					error = "difference raster size " + infiles[ifile] ;
					if (dataBuffer) delete[] dataBuffer;
					if (outBuffer) delete[] outBuffer;
					if (countBuffer) delete[] countBuffer;
					if (vcntbuffer) delete[] vcntbuffer;
					GDALClose(tds);
					
					return false ;
				}
			}
			else
			{
				GDALDataset* vcds = 0 ;
				if( useVcnt )
				{
					vcds = (GDALDataset*)GDALOpen( cntfiles[ifile].c_str(), GA_ReadOnly);
					if( vcds == 0 )
					{
						error = "Failed to open " + cntfiles[ifile] ;
						if (dataBuffer) delete[] dataBuffer;
						if (outBuffer) delete[] outBuffer;
						if (countBuffer) delete[] countBuffer;
						if (vcntbuffer) delete[] vcntbuffer;
						GDALClose(tds);
						return false ;
					}
					vcds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
							vcntbuffer, xsize, ysize, GDT_Int16 , 0, 0, 0);
					GDALClose(vcds) ;
				}

				tds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
					dataBuffer, xsize, ysize, GDT_Float32, 0, 0, 0);

				if( useVcnt )
				{
					for (int it = 0; it < asize; ++it)
					{
						if( vcntbuffer[it] > 0 )
						{
							//bugfixed 2018-4-13
							if( countBuffer[it] == 0 )
							{
								outBuffer[it]    = vcntbuffer[it] * dataBuffer[it] ;
								countBuffer[it]  = vcntbuffer[it] ;
							}else
							{
								outBuffer[it] += vcntbuffer[it] * dataBuffer[it] ;
								countBuffer[it] += vcntbuffer[it] ;
							}
						}else if( countBuffer[it] == 0 )
						{
							outBuffer[it] = dataBuffer[it] ;
						}
					}
				}else
				{
					for (int it = 0; it < asize; ++it)
					{
						if (countBuffer[it] == 0)
						{
							outBuffer[it] = dataBuffer[it];
							if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
							{
								countBuffer[it]+=1;
							}
						}
						else
						{
							if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
							{
								outBuffer[it] += dataBuffer[it];
								countBuffer[it] += 1;
							}
						}
					}
				}
			}
			GDALClose(tds);
		}
		else
		{
			std::cout << "Failed to open " + infiles[ifile] << std::endl;
			if( ignorebadfile == false )
			{
				if (dataBuffer) delete[] dataBuffer;
				if (outBuffer) delete[] outBuffer;
				if (countBuffer) delete[] countBuffer;
				if (vcntbuffer) delete[] vcntbuffer;
				error = "failed to open " + infiles[ifile] ;
				return false ;
			}
		}
		std::cout << ifile << "/" << infiles.size() << std::endl;
	}

	for (int it = 0; it < asize; ++it)
	{
		if (countBuffer[it] > 0)
		{
			outBuffer[it] = outBuffer[it] / countBuffer[it];
		}
	}

	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (driver == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		if (vcntbuffer) delete[] vcntbuffer;
		error = "Failed to get driver of tiff.";
		return false;
	}
	GDALDataset* outds = driver->Create(outfile.c_str(), xsize, ysize, 1, datatype, 0);

	if (outds == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		if (vcntbuffer) delete[] vcntbuffer;
		error = "Failed to create out put : "+outfile ;
		return false;
	}

	outds->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,outBuffer , xsize,ysize,GDT_Float32 , 0,0,0);
	{
		outds->SetProjection(prjbuff) ;
		outds->SetGeoTransform(trans) ;
	}

	//记录有效个数
	{
		string cntfile = outfile + ".vcnt.tif" ;
		GDALDataset* outds1 = driver->Create(cntfile.c_str(), xsize, ysize, 1, GDT_Int16, 0);
		if (outds1 == 0)
		{
			if (dataBuffer) delete[] dataBuffer;
			if (outBuffer) delete[] outBuffer;
			if (countBuffer) delete[] countBuffer;
			if (vcntbuffer) delete[] vcntbuffer;
			error = "Failed to create out put : "+cntfile ;
			return false;
		}
		outds1->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,countBuffer , xsize,ysize,GDT_Int16 , 0,0,0);
		{
			outds1->SetProjection(prjbuff) ;
			outds1->SetGeoTransform(trans) ;
		}
		GDALClose(outds1) ;
	}
	
	delete[] dataBuffer; dataBuffer = 0;
	delete[] countBuffer; countBuffer = 0;
	delete[] outBuffer; outBuffer = 0;
	if (vcntbuffer) delete[] vcntbuffer;

	GDALClose(outds);
	std::cout << "combine ok." << std::endl;
	return true;
}



//2018-4-18 有效个数得以保存 指定一个vector用于说明对应files是否有对应vcnt，
//如果为false那么不检查这个文件是否有vcnt，继续计算
//2018-4-18
bool wgt_CoorConverter::datacombineAve3(
	std::vector<std::string>& infiles, std::string outfile, 
	double valid0, double valid1, std::string& error, bool ignorebadfile,
	vector<bool>& hasvcntVec)
{
	cout<<"wgt_CoorConverter::datacombineAve v3.0"<<endl ;
	std::cout << "starting combine for " << outfile <<" ignorebadfile:"<<ignorebadfile << std::endl;

	vector<string> cntfiles ;
	for( int i = 0 ; i<infiles.size() ; ++ i )
	{
		string cfile1 = infiles[i] + ".vcnt.tif" ;
		if( test_file_exists(cfile1) )
		{
			cntfiles.push_back(cfile1) ;
			cout<<infiles[i]<<" with vcnt "<<endl ;
		}else if( hasvcntVec[i] == false )
		{
			cntfiles.push_back("") ;
			cout<<infiles[i]<<" no vcnt "<<endl ;
		}else
		{
			cout<<infiles[i]<<" error no vcnt "<<endl ;
		}
	}

	bool useVcnt = false ;
	if( cntfiles.size() == 0 )
	{
		useVcnt = false ;
		cout<<"not using vcnt."<<endl ;
	}else 
	{
		if( cntfiles.size() == infiles.size() )
		{
			useVcnt = true ;
			cout<<"using vcnt."<<endl ;
		}else
		{
			error = "Error: not all input files have a .vcnt.tif file.";
			return false;
		}
	}


	//get output parameters.
	int xsize(0), ysize(0);
	GDALDataType datatype;
	double trans[6];
	char prjbuff[2048];
	for (int i = 0; i < infiles.size(); ++i)
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen(infiles[i].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			xsize = tds->GetRasterXSize();
			ysize = tds->GetRasterYSize();
			datatype = tds->GetRasterBand(1)->GetRasterDataType();
			strcpy(prjbuff, tds->GetProjectionRef());
			tds->GetGeoTransform(trans);
			GDALClose(tds);
			break;
		}
	}

	if (xsize == 0)
	{
		error = "Failed to get raster info,proj,geotrans.";
		return false;
	}

	int asize = xsize * ysize;
	float* dataBuffer = new float[asize];
	short* countBuffer = new short[asize];
	float* outBuffer = new float[asize];
	short* vcntbuffer = 0 ;
	if( useVcnt )
	{
		vcntbuffer = new short[asize];
		if( vcntbuffer==0 )
		{
			error = "Failed to alloc memories for vcntbuffer.";
			if (dataBuffer) delete[] dataBuffer;
			if (outBuffer) delete[] outBuffer;
			if (countBuffer) delete[] countBuffer;
			
			return false;
		}
	}

	if (dataBuffer == 0
		|| countBuffer == 0
		|| outBuffer == 0 )
	{
		error = "Failed to alloc memories.";
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		if (vcntbuffer) delete[] vcntbuffer;
		
		return false;
	}

	memset(countBuffer, 0, asize * sizeof(short));
	for (int ifile = 0; ifile < infiles.size() ; ++ifile )
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen( infiles[ifile].c_str(), GA_ReadOnly);

		if (tds != 0 )
		{
			int nx = tds->GetRasterXSize();
			int ny = tds->GetRasterYSize();
			if (nx != xsize || ny != ysize) {
				std::cout << "*** Error : " << infiles[ifile] << " raster data size is not equal with others. out." << std::endl;
				if( ignorebadfile == false )
				{
					error = "difference raster size " + infiles[ifile] ;
					if (dataBuffer) delete[] dataBuffer;
					if (outBuffer) delete[] outBuffer;
					if (countBuffer) delete[] countBuffer;
					if (vcntbuffer) delete[] vcntbuffer;
					GDALClose(tds);
					
					return false ;
				}
			}
			else
			{
				GDALDataset* vcds = 0 ;
				if( useVcnt )
				{
					if( hasvcntVec[ifile] )
					{
						vcds = (GDALDataset*)GDALOpen( cntfiles[ifile].c_str(), GA_ReadOnly);
						if( vcds == 0 )
						{
							error = "Failed to open " + cntfiles[ifile] ;
							if (dataBuffer) delete[] dataBuffer;
							if (outBuffer) delete[] outBuffer;
							if (countBuffer) delete[] countBuffer;
							if (vcntbuffer) delete[] vcntbuffer;
							GDALClose(tds);
							return false ;
						}
						vcds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
								vcntbuffer, xsize, ysize, GDT_Int16 , 0, 0, 0);
						GDALClose(vcds) ;
					}
				}

				tds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
					dataBuffer, xsize, ysize, GDT_Float32, 0, 0, 0);

				if( useVcnt && hasvcntVec[ifile] )
				{
					for (int it = 0; it < asize; ++it)
					{
						if( vcntbuffer[it] > 0 )
						{
							//bugfixed 2018-4-13
							if( countBuffer[it] == 0 )
							{
								outBuffer[it]    = vcntbuffer[it] * dataBuffer[it] ;
								countBuffer[it]  = vcntbuffer[it] ;
							}else
							{
								outBuffer[it] += vcntbuffer[it] * dataBuffer[it] ;
								countBuffer[it] += vcntbuffer[it] ;
							}
						}else if( countBuffer[it] == 0 )
						{
							outBuffer[it] = dataBuffer[it] ;
						}
					}
				}else
				{
					for (int it = 0; it < asize; ++it)
					{
						if (countBuffer[it] == 0)
						{
							outBuffer[it] = dataBuffer[it];
							if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
							{
								countBuffer[it]+=1;
							}
						}
						else
						{
							if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
							{
								outBuffer[it] += dataBuffer[it];
								countBuffer[it] += 1;
							}
						}
					}
				}
			}
			GDALClose(tds);
		}
		else
		{
			std::cout << "Failed to open " + infiles[ifile] << std::endl;
			if( ignorebadfile == false )
			{
				if (dataBuffer) delete[] dataBuffer;
				if (outBuffer) delete[] outBuffer;
				if (countBuffer) delete[] countBuffer;
				if (vcntbuffer) delete[] vcntbuffer;
				error = "failed to open " + infiles[ifile] ;
				return false ;
			}
		}
		std::cout << ifile << "/" << infiles.size() << std::endl;
	}

	for (int it = 0; it < asize; ++it)
	{
		if (countBuffer[it] > 0)
		{
			outBuffer[it] = outBuffer[it] / countBuffer[it];
		}
	}

	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (driver == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		if (vcntbuffer) delete[] vcntbuffer;
		error = "Failed to get driver of tiff.";
		return false;
	}
	GDALDataset* outds = driver->Create(outfile.c_str(), xsize, ysize, 1, datatype, 0);

	if (outds == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		if (vcntbuffer) delete[] vcntbuffer;
		error = "Failed to create out put : "+outfile ;
		return false;
	}

	outds->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,outBuffer , xsize,ysize,GDT_Float32 , 0,0,0);
	{
		outds->SetProjection(prjbuff) ;
		outds->SetGeoTransform(trans) ;
	}

	//记录有效个数
	{
		string cntfile = outfile + ".vcnt.tif" ;
		GDALDataset* outds1 = driver->Create(cntfile.c_str(), xsize, ysize, 1, GDT_Int16, 0);
		if (outds1 == 0)
		{
			if (dataBuffer) delete[] dataBuffer;
			if (outBuffer) delete[] outBuffer;
			if (countBuffer) delete[] countBuffer;
			if (vcntbuffer) delete[] vcntbuffer;
			error = "Failed to create out put : "+cntfile ;
			return false;
		}
		outds1->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,countBuffer , xsize,ysize,GDT_Int16 , 0,0,0);
		{
			outds1->SetProjection(prjbuff) ;
			outds1->SetGeoTransform(trans) ;
		}
		GDALClose(outds1) ;
	}
	
	delete[] dataBuffer; dataBuffer = 0;
	delete[] countBuffer; countBuffer = 0;
	delete[] outBuffer; outBuffer = 0;
	if (vcntbuffer) delete[] vcntbuffer;

	GDALClose(outds);
	std::cout << "combine ok." << std::endl;
	return true;
}


//2018-1-22
bool wgt_CoorConverter::datacombineMax(std::vector<std::string>& infiles, 
	std::string outfile, double valid0, double valid1, std::string& error , bool ignorebadfile )
{
	std::cout << "starting max combine for " << outfile <<" ignorebadfile:"<<ignorebadfile << std::endl;
	//get output parameters.
	int xsize(0), ysize(0);
	GDALDataType datatype;
	double trans[6];
	char prjbuff[2048];
	for (int i = 0; i < infiles.size(); ++i)
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen(infiles[i].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			xsize = tds->GetRasterXSize();
			ysize = tds->GetRasterYSize();
			datatype = tds->GetRasterBand(1)->GetRasterDataType();
			strcpy(prjbuff, tds->GetProjectionRef());
			tds->GetGeoTransform(trans);
			GDALClose(tds);
			break;
		}
	}

	if (xsize == 0)
	{
		error = "Failed to get raster info,proj,geotrans.";
		return false;
	}

	int asize = xsize * ysize;
	float* dataBuffer = new float[asize];
	short* countBuffer = new short[asize];
	float* outBuffer = new float[asize];
	if (dataBuffer == 0
		|| countBuffer == 0
		|| outBuffer == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to alloc memories.";
		return false;
	}

	memset(countBuffer, 0, asize * sizeof(short));
	for (int ifile = 0; ifile < infiles.size() ; ++ifile )
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen( infiles[ifile].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			int nx = tds->GetRasterXSize();
			int ny = tds->GetRasterYSize();
			if (nx != xsize || ny != ysize) {
				std::cout << "*** Error : " << infiles[ifile] << " raster data size is not equal with others. out." << std::endl;
				if( ignorebadfile == false )
				{
					if (dataBuffer) delete[] dataBuffer;
					if (outBuffer) delete[] outBuffer;
					if (countBuffer) delete[] countBuffer;
					GDALClose(tds);
					error = "difference raster size " + infiles[ifile] ;
					return false ;
				}
			}
			else
			{
				tds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
					dataBuffer, xsize, ysize, GDT_Float32, 0, 0, 0);
				for (int it = 0; it < asize; ++it)
				{
					if (countBuffer[it] == 0)
					{
						outBuffer[it] = dataBuffer[it];
						if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
						{
							countBuffer[it]+=1;
						}
					}
					else
					{
						if (dataBuffer[it] >= valid0 && dataBuffer[it] <= valid1)
						{
							if( outBuffer[it] < dataBuffer[it] )
							{
								outBuffer[it]=dataBuffer[it] ;
							}
						}
					}
				}
			}
			GDALClose(tds);
		}
		else
		{
			std::cout << "Failed to open " + infiles[ifile] << std::endl;
			if( ignorebadfile == false )
			{
				if (dataBuffer) delete[] dataBuffer;
				if (outBuffer) delete[] outBuffer;
				if (countBuffer) delete[] countBuffer;
				error = "failed to open " + infiles[ifile] ;
				return false ;
			}
		}
		std::cout << ifile << "/" << infiles.size() << std::endl;
	}

	for (int it = 0; it < asize; ++it)
	{
		if (countBuffer[it] > 0)
		{
			outBuffer[it] = outBuffer[it] / countBuffer[it];
		}
	}

	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (driver == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to get driver of tiff.";
		return false;
	}
	GDALDataset* outds = driver->Create(outfile.c_str(), xsize, ysize, 1, datatype, 0);
	if (outds == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (outBuffer) delete[] outBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to create out put : "+outfile ;
		return false;
	}

	outds->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,outBuffer , xsize,ysize,GDT_Float32 , 0,0,0);
	{
		outds->SetProjection(prjbuff) ;
		outds->SetGeoTransform(trans) ;
	}

	delete[] dataBuffer; dataBuffer = 0;
	delete[] countBuffer; countBuffer = 0;
	delete[] outBuffer; outBuffer = 0;

	GDALClose(outds);
	std::cout << "combine ok." << std::endl;
	return true;
}




//2018-4-11 有效数据个数求和
bool wgt_CoorConverter::datacombineSum(std::vector<std::string>& infiles, 
	std::string outfile, std::string& error , bool ignorebadfile )
{
	std::cout << "starting sum combine for " << outfile <<" ignorebadfile:"<<ignorebadfile << std::endl;
	//get output parameters.
	int xsize(0), ysize(0);
	GDALDataType datatype;
	double trans[6];
	char prjbuff[2048];
	for (int i = 0; i < infiles.size(); ++i)
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen(infiles[i].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			xsize = tds->GetRasterXSize();
			ysize = tds->GetRasterYSize();
			datatype = tds->GetRasterBand(1)->GetRasterDataType();
			strcpy(prjbuff, tds->GetProjectionRef());
			tds->GetGeoTransform(trans);
			GDALClose(tds);
			break;
		}
	}

	if (xsize == 0)
	{
		error = "Failed to get raster info,proj,geotrans.";
		return false;
	}

	int asize = xsize * ysize;
	short* dataBuffer = new short[asize];
	short* countBuffer = new short[asize];
	if (dataBuffer == 0
		|| countBuffer == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to alloc memories.";
		return false;
	}

	memset(countBuffer, 0, asize * sizeof(short));
	for (int ifile = 0; ifile < infiles.size() ; ++ifile )
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen( infiles[ifile].c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			int nx = tds->GetRasterXSize();
			int ny = tds->GetRasterYSize();
			if (nx != xsize || ny != ysize) {
				std::cout << "*** Error : " << infiles[ifile] << " raster data size is not equal with others. out." << std::endl;
				if( ignorebadfile == false )
				{
					if (dataBuffer) delete[] dataBuffer;
					if (countBuffer) delete[] countBuffer;
					GDALClose(tds);
					error = "difference raster size " + infiles[ifile] ;
					return false ;
				}
			}
			else
			{
				tds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
					dataBuffer, xsize, ysize, GDT_Int16, 0, 0, 0);
				for (int it = 0; it < asize; ++it)
				{
					countBuffer[it] += dataBuffer[it] ;
				}
			}
			GDALClose(tds);
		}
		else
		{
			std::cout << "Failed to open " + infiles[ifile] << std::endl;
			if( ignorebadfile == false )
			{
				if (dataBuffer) delete[] dataBuffer;
				if (countBuffer) delete[] countBuffer;
				error = "failed to open " + infiles[ifile] ;
				return false ;
			}
		}
		std::cout << ifile << "/" << infiles.size() << "\t" ;
	}
	std::cout<<std::endl ;

	GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
	if (driver == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to get driver of tiff.";
		return false;
	}
	GDALDataset* outds = driver->Create(outfile.c_str(), xsize, ysize, 1, GDT_Int16, 0);
	if (outds == 0)
	{
		if (dataBuffer) delete[] dataBuffer;
		if (countBuffer) delete[] countBuffer;
		error = "Failed to create out put : "+outfile ;
		return false;
	}

	outds->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,countBuffer , xsize,ysize,GDT_Int16 , 0,0,0);
	{
		outds->SetProjection(prjbuff) ;
		outds->SetGeoTransform(trans) ;
	}

	delete[] dataBuffer; dataBuffer = 0;
	delete[] countBuffer; countBuffer = 0;

	GDALClose(outds);
	std::cout << "combine ok." << std::endl;
	return true;
}




//2018-1-28
ccRasterInfo wgt_CoorConverter::getRasterInfo(std::string filepath)
{
	ccRasterInfo ri;
	ri.xsize = 0;
	ri.ysize = 0;
	GDALDataset* ds = (GDALDataset*)GDALOpen(filepath.c_str(), GA_ReadOnly);
	if (ds != 0)
	{
		ri.xsize = ds->GetRasterXSize();
		ri.ysize = ds->GetRasterYSize();
		ds->GetGeoTransform(ri.trans);

		OGRSpatialReference osr;
		char prjbuff[2048];
		char* ptrPrj = prjbuff;
		char* ptr2 = 0;
		strcpy(prjbuff, ds->GetProjectionRef());
		osr.importFromWkt(&ptrPrj);
		osr.exportToProj4(&ptr2);

		ri.projstr = std::string(prjbuff);
		ri.proj4str = std::string(ptr2);

		CPLFree(ptr2); ptr2 = 0;

		GDALClose(ds);
	}
	return ri;
}



//2018-2-12
bool wgt_CoorConverter::dataSubtract( std::string fileA , std::string fileB , std::string outfile , 
	ccDataInfo ia, ccDataInfo ib, double outfillvalue, std::string& error )
{
	std::cout << "starting subtract : " << fileA << " - " << fileB << std::endl;
	std::cout << "result : " << outfile << std::endl;
	//get output parameters.
	int xsize(0), ysize(0);
	double trans[6];
	char prjbuff[2048];
	{
		GDALDataset* tds = (GDALDataset*)GDALOpen(fileA.c_str(), GA_ReadOnly);
		if (tds != 0)
		{
			xsize = tds->GetRasterXSize();
			ysize = tds->GetRasterYSize();
			strcpy(prjbuff, tds->GetProjectionRef());
			tds->GetGeoTransform(trans);
			GDALDataset* bds = (GDALDataset*)GDALOpen(fileB.c_str(), GA_ReadOnly);
			if( bds != 0 )
			{
				int bxsize = bds->GetRasterXSize() ;
				int bysize = bds->GetRasterYSize() ;
				if( xsize == bxsize && ysize == bysize )
				{
					int asize = xsize * ysize;
					float* bufferA = new float[asize];
					float* bufferB = new float[asize];
					if (bufferA == 0
						|| bufferB == 0)
					{
						GDALClose(tds);
						GDALClose(bds) ;
						if (bufferA) delete[] bufferA;
						if (bufferB) delete[] bufferB;
						error = "Failed to alloc memories.";
						return false;
					}else
					{
						tds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
							bufferA, xsize, ysize, GDT_Float32, 0, 0, 0);
						bds->GetRasterBand(1)->RasterIO(GF_Read, 0, 0, xsize, ysize ,
							bufferB, xsize, ysize, GDT_Float32, 0, 0, 0);
						GDALClose(tds);
						GDALClose(bds) ;
						for (int it = 0; it < asize; ++it)
						{
							if (
								bufferA[it] >= ia.valid0 && bufferA[it] <= ia.valid1
								&& bufferB[it]>=ib.valid0 && bufferB[it]<=ib.valid1 
								)
							{
								bufferA[it] = (bufferA[it]*ia.scale+ia.offset)-(bufferB[it]*ib.scale+ib.offset);
							}else
							{
								bufferA[it] = outfillvalue ;
							}
						}
					}
					GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
					GDALDataset* outds = driver->Create(outfile.c_str(), xsize, ysize, 1, GDT_Float32, 0);
					if (outds == 0)
					{
						if (bufferA) delete[] bufferA;
						if (bufferB) delete[] bufferB;
						error = "Failed to create output : "+outfile ;
						return false;
					}
					outds->SetProjection(prjbuff) ;
					outds->SetGeoTransform(trans) ;
					outds->GetRasterBand(1)->RasterIO(GF_Write,0,0,xsize,ysize,bufferA , xsize,ysize,GDT_Float32 , 0,0,0);
					GDALClose(outds);
					if (bufferA) delete[] bufferA;
					if (bufferB) delete[] bufferB;
					std::cout << "done." << std::endl;
					return true;
				}else
				{
					GDALClose(tds);
					GDALClose(bds) ;
					error = std::string("DataA size is different with DataB") ; 
					return false ;
				}
			}else
			{
				GDALClose(tds);
				error = std::string("Failed to open ") + fileB ; 
				return false ;
			}
		}else
		{
			error = std::string("Failed to open ") + fileA ; 
			return false ;
		}
	}
}


#endif
